bulk-org
Raphael Kubler
7/17/2023
Introduction
Bulk RNA-seq analysis of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).
CO vs COiMG
Organoids co-cultured with vs without microglia.
Variance analysis
PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
theme_classic()
Heatmap
#for control only
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'
#horizontal
ra <- rowAnnotation(
COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
col = list(
COiMG = c("COiMg"='tomato','CO'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr <- t(batch.rem3_mod1[,meta3$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
DEA
Volcano
volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
res <-
mutate(res,
sig = case_when(
padj >= 0.05 ~ "non_sig",
padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
)) %>%
mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
mutate(log2FoldChange = case_when(
log2FoldChange > 10 ~ Inf,
log2FoldChange < -10 ~ -Inf,
TRUE ~ log2FoldChange
)) %>%
mutate(class = paste(sig, direction))
if( type == "ALS"){
xpos <- 0.5
ymax <- ymax1
xlim <- c(xmax1,xmax2)
}else{
xpos <- 0.025
ymax <- ymax2
xlim <- c(-0.042, 0.042)
}
de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
filter(sig != "non_sig") %>%
mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
mutate(n = formatC(n, format="f", big.mark=",", digits=0))
plot <- res %>%
mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) +
#geom_point(aes(colour = class ), size = 0.5) +
rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
scale_colour_manual(values = c("non_sig up" = "gray",
"non_sig down" = "gray",
"sig up" = "#EB7F56",
"sig - strong up" = "#B61927",
"sig down" = "#4F8FC4",
"sig - strong down" = "dodgerblue4"
)) +
theme_bw() +
labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
guides(colour = FALSE) +
scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.border = element_blank(),
axis.ticks = element_line(colour = "black")
) +
geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
scale_x_continuous(limits = xlim)
if(!is.null(annotate_by)){
plot <- plot +
ggrepel::geom_text_repel(
fontface = "italic",
data = filter(res, symbol %in% annotate_by),
aes(x = log2FoldChange, y = -log10(pvalue), label = symbol),
min.segment.length = unit(0, "lines"),
size = 2.3) +
geom_point(
data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
) +
geom_point(aes(colour = class ),
data = filter(res, symbol %in% annotate_by), size = 0.6
)
}
return(plot)
}
res_COiMG$symbol <- rownames(res_COiMG)
volcano_plot(data.frame(res_COiMG), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top15.coimg <- c(rownames(results$COiMG[order(results$COiMG$log2FoldChange),][results$COiMG[order(results$COiMG$log2FoldChange),]$padj < 0.05,][1:25,]),
rownames(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,]))
ha <- HeatmapAnnotation(
COiMG = meta3$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')))
ComplexHeatmap::Heatmap(t(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top15.coimg,]))),
cluster_rows = T,
cluster_columns = T,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
show_column_names = F,
top_annotation = ha,
row_names_gp = gpar(fontsize=8),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(30),
name = "Z-score expression",
)
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
knitr::kable(cbind('DEG'=top15.coimg,'Cilia-associated'=top15.coimg %in% organoid$CBC))| DEG | Cilia-associated |
|---|---|
| CD164L2 | FALSE |
| STOML3 | TRUE |
| MMP10 | FALSE |
| GIPC2 | FALSE |
| LRRC18 | FALSE |
| FAM166C | FALSE |
| LRRC71 | FALSE |
| IL5RA | FALSE |
| TTC22 | FALSE |
| CD302 | FALSE |
| C20orf85 | FALSE |
| CLDN2 | FALSE |
| MFRP | FALSE |
| TACR3 | FALSE |
| CIBAR2 | FALSE |
| C11orf97 | FALSE |
| F5 | FALSE |
| SAG | FALSE |
| ADGB | FALSE |
| TEKT4 | FALSE |
| VWA5B1 | FALSE |
| DEUP1 | FALSE |
| DYNLT5 | FALSE |
| ARSI | FALSE |
| FAM81B | TRUE |
| APOC2 | FALSE |
| GOLGA6D | FALSE |
| VAV1 | FALSE |
| FCGR3A | FALSE |
| C1QC | FALSE |
| ITGB2 | FALSE |
| CD86 | FALSE |
| OLR1 | FALSE |
| SPP1 | FALSE |
| WDFY4 | FALSE |
| ALOX5AP | FALSE |
| TYROBP | FALSE |
| ADAMDEC1 | FALSE |
| BIN2 | FALSE |
| CD53 | FALSE |
| MS4A7 | FALSE |
| SRGN | FALSE |
| CCRL2 | FALSE |
| CD52 | FALSE |
| FCER1G | FALSE |
| LSP1 | FALSE |
| PECAM1 | FALSE |
| CD48 | FALSE |
| DCSTAMP | FALSE |
| CHI3L1 | FALSE |
GO
dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$COiMG$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts3)]), check.names = F)
for(i in names(cgenes)[-1]){
k <- na.omit(cgenes[i])
s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts3),]
m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
IFNy stimulation
Let’s do this for IFNy stimulation
Variance analysis
PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
#stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
theme_classic()
DEA
Volcano
res_IFNg$symbol <- rownames(res_IFNg)
volcano_plot(data.frame(res_IFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
COiMg.degs <- all_res$COiMG$DEGsTop 50 DEGs heatmap
top15.ifng <- c(rownames(results$IFNg[order(results$IFNg$log2FoldChange),][results$IFNg[order(results$IFNg$log2FoldChange),]$padj < 0.05,][1:25,]),
rownames(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,]))
ha <- HeatmapAnnotation(Condition = meta3$condition,
COiMG = meta3$COiMg,
Line = meta3$Line,
col = list( Condition = c("ctrl" = "lightblue", "IFNg" = "tomato", "LPS" = "darkred"),
COiMG = c("yes"='darkblue','no'='grey80'),
Line = c('WTC11'='black','MSN38'='white')))
ComplexHeatmap::Heatmap(t(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top15.ifng,]))),
cluster_rows = T,
cluster_columns = T,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
show_column_names = F,
top_annotation = ha,
row_names_gp = gpar(fontsize=8),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(30),
name = "Z-score expression",
)
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
knitr::kable(cbind('DEG'=top15.ifng,'Cilia-associated'=c(top15.ifng %in% organoid$CBC)))| DEG | Cilia-associated |
|---|---|
| CROCC2 | FALSE |
| F5 | FALSE |
| STMND1 | FALSE |
| MUSK | FALSE |
| ARHGAP36 | FALSE |
| PTPRC | TRUE |
| IGF1 | FALSE |
| KL | FALSE |
| STOML3 | TRUE |
| ATP13A5 | FALSE |
| SLC39A12 | TRUE |
| CYTL1 | FALSE |
| C11orf97 | FALSE |
| SERPINA3 | FALSE |
| TMEM125 | FALSE |
| FAM81B | TRUE |
| C20orf85 | FALSE |
| CAPSL | TRUE |
| CIBAR2 | FALSE |
| TRPM1 | FALSE |
| C10orf105 | FALSE |
| HMGCS2 | FALSE |
| CFAP73 | FALSE |
| PMCH | FALSE |
| MMP10 | FALSE |
| IDO1 | FALSE |
| GBP5 | FALSE |
| ETV7 | FALSE |
| CXCL10 | FALSE |
| SIGLEC7 | FALSE |
| GBP1 | FALSE |
| CXCL9 | FALSE |
| APOL6 | FALSE |
| APOL1 | FALSE |
| PLAAT4 | FALSE |
| GBP3 | FALSE |
| SAMD9L | FALSE |
| GBP2 | FALSE |
| XAF1 | FALSE |
| PSMB9 | FALSE |
| CXCL11 | FALSE |
| GBP6 | FALSE |
| CD274 | FALSE |
| OAS2 | FALSE |
| DES | FALSE |
| UBD | FALSE |
| OR2I1P | FALSE |
| OAS1 | FALSE |
| NLRC5 | FALSE |
| BATF2 | FALSE |
GO
dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
For organoid cell-types.
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For NDD gene lists.
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]## [1] "CYP7B1" "GRIN2A"For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)